clear all; close all; clc
 addpath('../LibForSerial/LibrarySLP','../LibForSerial/LibraryLeSage/');
%----------- Load shocks using two identification schemes ----------------%
  load('./DataShocks/dataJMNupdated.mat');  
 uF        = ehat_maxG(:,3); clearvars -except uF Uf Um datem 

%---------- Construct indicators for seq. of shocks with LMN ident scheme-%
IndPreviousNShocks       = zeros(size(uF));
  NumberConsecutiveshocks = 2; % 1 previous shocks that is  positive, and the 2nd one is positive too 
%/////////////////////////////////////////////////////////////////////////
for tt=NumberConsecutiveshocks:length(uF)-1 
    if tt==NumberConsecutiveshocks
        if (                                    sum(uF(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks                  );  IndPreviousNShocks(tt) = 1; end
    else
        if (uF(tt-NumberConsecutiveshocks)<0 && sum(uF(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks                  );  IndPreviousNShocks(tt) = 1; end
    end
end

clear tt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sample=[ 1961 1; 2019  12];
data = xlsread('../Dataset/VARdataBBupdatedMonthly.xlsx','MainData');
startsample=find(data(:,1)==sample(1,1)*100+sample(1,2));
endsample  =find(data(:,1)==sample(2,1)*100+sample(2,2));

data = data(startsample:endsample,2:end);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DATA TRANSFORMATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data(:,[3:10]) = log(data(:,[3:10]))*100;
%--------------------------------------%
clear startsample endsample 

T = size(data,1);
%P  =  3;
 P  =  6; % including in the LP regression P lags of all variables in the system "use p = 6 lags in the VARs" as in LMN AEJ:Macro
%P  = 12; %                                                      " ... noting that using 12 lags makes no difference to the results."
 
  plottype = 'Output';      ii=3; %
%% IRF 
h1 = 0 ;  %start LP at h=0 or h=1 (h=1 if impose no contemporaneous impact)
H  = 12*3; %#horizon of IR

% We are interested in the estimation of the dynamic multiplier of y_t+h with respect to a change in x_t for h ranging from 0 to H
%------------ IRF of GDP growth to a shock variable ----------------------%
    y  = data(:,ii);    %endogenous variable
 if ii == 8 || ii == 9 
     y  = [NaN(12,1); data(13:end,ii)-data(1:end-12,ii)];    %endogenous variable
 end
%-------------------------------------------------------------------------%
Remove_Large_Shocks = 0;
    x  = uF;      
    IndSeqShocks = IndPreviousNShocks;        % JLN Financial Uncertainty shocks (from trivariate system with Macro Uncertainty)
    IndLargeShocks = uF>2;

LargeShocksToBeRemoved = IndSeqShocks & IndLargeShocks;  
% [x IndSeqShocks IndLargeShocks]
  if Remove_Large_Shocks; IndSeqShocks = IndSeqShocks - LargeShocksToBeRemoved; end
fprintf('Number of events with %d consecutive shocks for identification scheme with uM is %d \n',NumberConsecutiveshocks,sum(IndSeqShocks));
clear uF2015 IndPreviousNShocks IndPreviousNShocks2018 vardata

%/////////////////////////////////////////////////////////////////////////%
%-------------------------------------------------------------------------%
xs = x.*IndSeqShocks; %shock variable times states
%-------------------------------------------------------------------------%
 %w  = [lagmatrix( data(:,unique([ii 1 10])), 1:P )                       ]; %control variables (contemporaneous vars, lagged vars)
 %w  = [lagmatrix( data(:,unique([ii 1])   ), 1:P )                       ]; %remove employment since not in the first step of the LMN AEJ Macro VAR
  w  = [lagmatrix( data(:,unique([ii 1])   ), 1:P ) lagmatrix(Uf,1:P) lagmatrix(Um,1:P)]; %Replicating LMN AEJ:Macro VAR system (removing the RF changes little)

%/////////////////////////////////////////////////////////////////////////%  
% Additional robustness not reported in the paper; RUN with benchmark case P=6; 
%/////////////////////////////////////////////////////////////////////////%
%%w  = [lagmatrix( data(:,unique([ii 1 6])), 1:P )                       ];
%%w  = [lagmatrix( data(:,unique([ii 1  ])), 1:P )       data(:,6)       ];
 
  w( ~isfinite(w) ) = 0;
OKforReg = isfinite(x) & isfinite(y) & isfinite(xs);
for cc=1:size(w,2)
    OKforReg = OKforReg & isfinite(w(:,cc)) ;
end  
  x = x(OKforReg,:); 
  y = y(OKforReg,:); 
  w = w(OKforReg,:);
  xs=xs(OKforReg,:);
  
% This shrinks the estimated dynamic multiplier to a polynomial of order r - 1
r = 3; %(r-1)=order of the limit polynomial (so r=3 implies the IR is shrunk towards a quadratic polynomial)
lambda = 20; %value of penalty
if ii==8 || ii==9;        lambda = 10000; end                                      
if ii==6 || ii==7; r = 2; lambda = 10000; end  % r=2 mostly to smooth the oscillations at h=36 months, for which there is huge uncertainty anyway!                 
%-------------------------------------------------------------------------%
% k-fold cross validation 
  k_fold=5;
%-------------------------------------------------------------------------%
slp_sd  = locprojStateDep(y,x                  ,xs,w,h1,H,'smooth',r,lambda); %IR from Smooth Local Projection
 lp_sd  = locprojStateDep(y,x                  ,xs,w,h1,H,'reg');             %IR from        Local Projection

%% Figure 4, Panel (a) State-dependent IR to a FIN shock
figure(4)
%-------------------------------------------------------------------------%
     fname = sprintf('./DataShocks/LudvigsonMaNg_LP_linear_0.mat');
load(fname)
if       ii==1;             slp_linear = slpShortRate;   slp_linearUp = slpShortRateUP;   slp_linearDown = slpShortRateDown; 
elseif   ii==3;             slp_linear = slpOutput;      slp_linearUp = slpOutputUP;      slp_linearDown = slpOutputDown;
elseif   ii==6 || ii==7;    slp_linear = slpStockMkt;    slp_linearUp = slpStockMktUP;    slp_linearDown = slpStockMktDown;
elseif   ii==5 || ii==9;    slp_linear = slpInflation;   slp_linearUp = slpInflationUP;   slp_linearDown = slpInflationDown;
end
%-------------------------------------------------------------------------%
    sigma_state = 1;
if   ii==9 
    sigma_state = 0.5;
end

%-------------------------------------------------------------------------%
subplot(1,2,1)
if   ii==1 || ii==3 || ii==7 || ii==9
  plot( 0:H ,slp_linear, 'k', 0:H, (slp_sd.IR+slp_sd.IRstate*sigma_state), 'b--', 0:H, (slp_sd.IR),'k--', 'LineWidth', 1.5 ); 
else
  plot(                       0:H, (slp_sd.IR+slp_sd.IRstate*sigma_state), 'b--', 0:H, (slp_sd.IR),'k--', 'LineWidth', 1.5 ); 
end
  title(['IR_{slp} when there are ' num2str(NumberConsecutiveshocks) ' Consecutive Positive Shocks']); %legend(          's_t=+1',                      's_t=0','Location','Best')
hold on; plot(0:H , zeros(H+1,1) , '-k' , 'LineWidth' , 2 ); grid; xlim([0 H]);                          legend('Linear','s_t=+1 (Previous shock(s)=1)','s_t=0','Location','Best')
xlim([0      H]);set(gca,'XTick',[0 3:3:H],'FontSize',12);

subplot(1,2,2)
grpyat=[(0:1:H)', slp_sd.Interval_up_state*sigma_state; (H:-1:0)' slp_sd.Interval_down_state(end:-1:1)*sigma_state]; patch(grpyat(:,1), grpyat(:,2), [0.7 0.7 0.7],'edgecolor', [0.7 0.7 0.7]); hold on
plot(0  :H , zeros(H+1,1), 'k-'); hold on 
plot(    0:1:H,   slp_sd.IRstate*sigma_state, 'b--','LineWidth', 1.5); hold on
title('State multiplier of State-dependent IR to uncertainty shock')
grid;
xlim([0      H]);set(gca,'XTick',[0 3:3:H],'FontSize',12);